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Abstract. We show that the phenomenon of modulational instability in arrays of 
Bose-Einstein condensates confined to optical lattices gives rise to coherent spatial 
structures of localized excitations. These excitations represent thin disks in ID, narrow 
tubes in 2D, and small hollows in 3D arrays, filled in with condensed atoms of much 
greater density compared to surrounding array sites. Aspects of the developed pattern 
depend on the initial distribution function of the condensate over the optical lattice, 
corresponding to particular points of the Brillouin zone. The long-time behavior of 
the spatial structures emerging due to modulational instability is characterized by the 
periodic recurrence to the initial low-density state in a finite optical lattice. We propose 
a simple way to retain the localized spatial structures with high atomic concentration, 
which may be of interest for applications. Theoretical model, based on the multiple 
scale expansion, describes the basic features of the phenomenon. Results of numerical 
simulations confirm the analytical predictions. 
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1. Introduction 

Optical lattices formed by laser waves are the media where Bose-Einstein condensates 
(BEC) exhibit remarkable properties. Over the last few years there has been a significant 
progress in manipulation with BEC confined to optical lattices which resulted in 
observation of diverse new phenomena, such as coherent emission of Bose-condensed 
atoms PP, Bloch oscillations and Landau-Zener tunneling |2], atomic Josephson effect j^j, 
Mott insulator - superfluid transition [I]. Realization of BEC arrays in 2D and 3D optical 
lattices jUEj has opened new perspectives for investigation of fundamental properties of 
quantum gases in lower dimensions. Theoretical studies on the dynamics of BEC in a 
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periodic potential critically rely on the concepts of band structure and Bloch states, 
acquired from the solid state physics. Recent developments in the field show that 
the above concepts, originally constructed for linear periodic systems, also play the 
crucial role in the physics of nonlinear periodic systems, such as BEC arrays in optical 
lattices jOdlE] • A good correlation between predictions of the band theory for the BEC 
dynamics in a periodic potential and experimental results with static, moving, and 
accelerating ID optical lattices was demonstrated in [0|. Wide range of spatiotemporal 
behaviour of a BEC in ID linear- and circular-chain optical lattices in the tight-binding 
limit, when the dynamics is described by the discrete nonlinear Schrodinger equation, 
was numerically investigated in Ref. [101 ■ Transmission of matter wave pulses incident 
in a ID optical lattice, including their collision dynamics, was theoretically considered 
indU- 

A subject, which is interesting both from the viewpoints of the theory and 
applications of BEC, is the origin and dynamics of spatially localized nonlinear 
excitations in the condensate confined to a periodic potential. A stimulating discovery 
was the proof of the existence of bright solitons in the effectively ID BEC arrays with 
repulsive interaction between atoms Pll2pi3j . In view of the fact that a continuous BEC 
with repulsive interatomic forces does not support spatially localized humps of atomic 
concentration, BEC arrays in optical lattices are considered to be the most inviting 
media for the creation and manipulation of soliton-like structures with ultracold atoms. 
The physical mechanism by which this possibility arises is similar to that of electrons in a 
periodic potential, in specific cases acquiring negative effective mass. The presence of the 
optical lattice can invert the sign of the dispersive term, which then balances the action 
of the nonlinearity. Therefore, bright solitons in BEC arrays with repulsive interaction 
between atoms are possible in the presence of the periodic potential of the optical 
lattice. It is worth to mention here that although a continuous BEC with attractive 
interaction between atoms can bear bright solitons, its other property leading to a 
collapse of the condensate at some critical atomic concentration (for review see e.g. |14j). 
makes it less appealing for the above purpose. In recent papers [El EE], reporting 
on the first experimental observation of matter-wave bright solitons in a continuous 
BEC with attractive interatomic forces ( 7 Li), the macroscopic quantum bound state 
of Bose-condensed atoms (bright soliton) was shown to exist in a narrow window of 
atomic numbers (around M ~ 5000). Beyond that window atomic wavepackets undergo 
collapse or explosion. Moreover, the bright soliton with attractive forces between atoms 
subject to expulsive potential, as applied in [TSjfTBj . appears to be of limited lifetime due 
to the effect of quantum tunneling [T2j (termed by authors as quantum evaporation) , 
which leads to eventual explosion of the soliton. Therefore, matter-wave bright solitons 
composed of repulsive atoms in optical lattices, which are free of the above constraints, 
seem to be advantageous for applications. 

Out of existing studies on localized excitations in BEC arrays, little attention has 
been devoted to methods of creation of such structures, so far. It has recently been 
suggested to employ the modulational instability, which constitutes one of the most 
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important phenomena associated with the evolution of nonlinear waves, to create bright 
BEC solitons in a ID optical lattice jH] . A variety of localized solutions are found to the 
one-dimensional nonlinear Schrodinger equation with a periodic potential, some of which 
are spatially and temporally stable [12]. Interesting consequence of the modulational 
instability in a continuous BEC was reported in jTH]. The authors have shown that the 
modulational instability leads to fragmentation of the ferromagnetic phase in a spinor 
Bose-condensate. Another manifestation of the modulational instability as leading to 
dynamical superfluid-insulator transition in a BEC confined to an optical lattice and 
magnetic potential has been studied in |T9"] . 

It is appropriate to mention, that the phenomenon of modulational instability is 
well studied in different areas of nonlinear physics, since initiated in the 1960s, by 
predictions in hydrodynamics [20J, plasmas [2lJ|22], and nonlinear optics [SHIEI]- For 
later reviews on modulational instability in Hamiltonian systems the reader is addressed 
to references (23123123 • I n view of the existence of many features of ultracold atomic 
gases similar to those observed in the above mentioned systems, there is a solid ground 
to expect rich dynamics induced by the modulational instability in such a nonlinear 
system as Bose-Einstein condensates. 

In the present paper we study the dynamical processes in BEC arrays confined to 
one-, two-, and three-dimensional optical lattices, which are due to the modulational 
instability Particularly we focus on the coherent spatial structures in 2D and 3D BEC 
arrays, which originate from the modulational instability. 

The paper is organized as follows: Section 121 contains the derivation of main 
equations, as well as brief exposition of the modulational instability. In section |3] we 
present the energy band structure for a BEC distributed over the periodic potential. In 
section 0] we analyze the features of spatial structures in arrays of BEC, originated form 
the modulational instability of the initial waveforms, corresponding to different points 
of the Brillouin zone. Section summarizes the results of this study. 



2. The multiscale analysis and modulational instability 

To develop the model we consider a dimensionless 3D Gross-Pitaevskii (GP) equation 

i^-A^ + y^ + xi^, (i) 

where r = (r x ,r y ,r z ). In (fl|) the spatial coordinates are normalized to £, I being a 
characteristic size of the potential (say, the smallest of its periods), the time is mesured 
in units of 2m£ 2 /h, and the amplitude of the order parameter is normalized to the total 
number of atoms per unit volume a/A/"/£ 3 . Then the nonlinearity coefficient \ is given by 
X = SnJ\fa s /£, where a s is the s-wave scattering length. The potential V(r) is assumed, 
for the sake of simplicity, to be separable, i.e. of the form V(r) = J2j ^j( r j)> 3 = x i Vi z 
(which corresponds to the majority of experimental settings), and periodic in each of 
the spatial directions: Vj(rj) = Vj(rj + aj), with aj the period in the direction rj 
(in accordence with the accepted scaling aj > 1). For convenience, the equation (JTJ) is 
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considered subject to periodic boundary conditions ^(r) = ip{r x + L X) r y , r z ), etc., where 
Lj = Njaj with Nj and Lj respectively, the number of primitive cells and the length of 
the system in the direction rj. The theory is developed for the small amplitude limit, 
when the multiscale analysis is applicable. Hence, we look for a solution to equation (0) 
in the form 

ip = eipi + e 2 ip 2 + e 3 ip 3 + ••• (2) 

where the ipj are functions of the scaled independent variables r p = e p t, £ p = e p r, 
p = 0, 1, 2, with e a small parameter to be specified later. Denoting with u aj (qj), and 
$ a .(r 3 ) = \oij,qj), the eigenvalues and eigenfunctions of the periodic operators C r . = 
— Qj? + Vj(rj), we have that the solution to a linear part of the equation ((H), Cip = 0, 
with £ = idt — ^2j£ rj , can be written in the form \m x m y m z ) = Y[j ■,( r i) eU ^°' i 
with <& m (rj) Bloch states of the corresponding ID linear problems. Here rrtj denotes the 
couple of quantum numbers {azj,qj}, with aj the band index and qj the component of 
the wave vector in the j direction (note that the imposed boundary conditions obviously 
imply that qj = qj^ n = j^n so that the extension of the Brillouin zone (BZ) in the j 
direction is [— 7r/aj, 7r/aj]). 

Substituting the equation (J2J) into (JTJ), and collecting terms of equal powers in e, 
one arrives at the set of equations Cip n = M n , with 



Mi = 0, M 2 = -%d Tx i\) X - 2V ViVi, 

M 3 = -id n ^ x - %3 Tx i) 2 - A^i - 2V • V X V>2 - 2V ■ V 2 ^i + xl^ilVi, 

where V p denotes the gradient with respect to £ p . 

Since we are interested in instabilities of the condensate wavefunction, we 
investigate the influence of the nonlinear term in the equation ([TJ on the Bloch states 
of the underlying linear problem. To this end we take as starting point in the expansion 
(J2J), a modulated state of the form 

V>i = A{ix, n, ...)e~ lu, ° T0 \m 0x m 0y m 0z ), (3) 

with co> = J2j UaojiQj) (the subscript zero refers to the chosen band, below we consider 
the two lowest ones). Then the first order equation is automatically satisfied by if)i, 
while the equation of the second order can be solved in the form 

ip 2 = ^2'B a e~ tuJ0T0 \a x , q 0>x ; a y , q 0tV ; a z , q 0>z ), (4) 

a 

where the prime denotes a ^ ao in the sum and have taken into account that the terms 
with q 7^ qo give zero contribution. Analysis similar to that of reference [H] shows that 
A = A(R; ^2--5 T 2, •••) with R = ^ — vr x and v = — (a 0x a 0y a 0z \2i\7\a 0x a 0y ao z ) is the 
group velocity of the carrier wave. The coefficients B a are found as 

r (,V,z) ft a, p(x,«) ft A + -V&'V) ft A 
Oa — / v > V>) 

uj - uj a {q ) 



Regular spatial structures in arrays of BEC 5 

where u; a (q ) = Z)j w «j(9bi)> r &ao,« = tfo^SJao,*, gb,*)*a I) ,ao, v *o„ao,. (other 

coefficients T are obtained by cyclic permutations of x,y,z). Finally, considering the 
orthogonality of to \mo x mo y mo z ) we obtain the following 3D NLS for the slowly 
varying envelope 

j=x,y,z 3 

where we assumed A not depending on £2, and introduced the inverse of the effective 
mass tensor 

and the effective nonlinearity 

x=x n / i $ m0J i 4 ^- (s) 

j=x,y,z 

Now we are at the point to discuss the small parameter. To simplify, we consider 
a cubic box with aj = I and Lj = L. On the one hand, as it was mentined above, 
the physical order parameter is normalized to the total number of atoms A/", while the 
formal wave function ip must be normalized to one: J \ip\ 2 dr = 1. On the other hand 
all parameters in the equation ((Hj) must be of order one. Next we notice the oscillatory 
character of the Bloch functions, in a general situation leads to the fact that the integrals 
in the last expression (jSJ) for x has a numerical smallness (see e.g. the examples below). 
Then, taking into account that x — SnAfa s /£ one can define e = Ma s £ 2 /L 3 . Consider 
now a condensate with M = 10 5 of 87 Rb atoms (a s ~ 5.5 nm) homogeneously distributed 
over a cubic box with L = 100 /im having iVj = 100 cells in each direction (respectively 
I = 1 /im). Then we compute e ~ 0.014. A physical situation when e is not too small (say 
in experiments 5J it can be identified ase~ 0.257) the multiple scale expansion, strictly 
speaking, is not valid. For this reason below we employ the numerical simulations, 
which, however, clearly illustrate that the small amplitude limit gives remarkably good 
estimates for the characteristic scales of the problem and allows one to understand the 
symmetry of the developed patterns. 

Let us analyze the stability problem within the framework of equation (jSJ, i.e. look 
for a solution of the form © with A = (p + ae i(nT2 ~ KR) + fr e -*(™2-KR)j e -Vr 2; w h ere 
|a| , |fe| <C p- This solution is modulationally unstable if 

Z{Z + Axp 2 ) < 0, Z = J2 (9) 

j=x,y,z 

where the modulational instability is understood in the sence of a plain wave instability 
with respect to small modulations of its amplitude. Below in section 0] we consider the 
modulational instability of solutions to equation (JTJ) in the form of Bloch states induced 
by periodic small amplitude and long-wavelength perturbation. 
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3. Energy band structures 

In the previous section we assumed the knowledge of the Bloch states and energy 
bands of the underlying linear Schrodinger problem. For generic multidimensional 
potentials this could be a quite difficult problem to solve. To avoid difficulties we 
shall restrict here to the case of separable potentials of trigonometric form i.e. we take 

^( r ) = Ej V j( r j), J = x ,V,z with V j( r j) g iven b Y 

Vj(rj) = Acos^r,). (10) 

Here A is a constant fixing the depth of the lattice and 2n/Kj the periodicity in the 
rj directions. In the following we fix Kj = 2 for all j so that the potential will be a 
superposition of identical one dimensional Mathieu potentials, and the corresponding 
Bravais lattice will have simple cubic symmetry. The band structure and the Bloch 
states of the linear Schrodinger equation are then obtained as 

S a (k) = ^e a (k,) l * Q (k,r) = n^(k j ,r j ), (11) 
j j 

where a denotes the band index and e a (kj) and i/; a (kj,rj) are the eigenvalues and 
eigenfunctions of the Mathieu equation 

T^T + A cos(2x)v9 fc>Q = E a (k)ip k)Cc . (12) 

This equation can be transformed, making use of the expansion of the wave function 
into momentum eigenfunctions (Fourier expansion), into a tridiagonal problem whose 
solutions can be obtained with high accuracy by means of continued fractions. As an 
example of these calculations we report in figure Q the first two energy bands of the 2D 
separable Mathieu potential, in the case A = 1 (note that with the choice k = 2 the 
potential has periodicity 7r in both directions, so that the BZ is a square of size 2). 

A=l 




Figure 1. The first two bands of the 2D separable Mathieu potential of amplitude 
A = 1 and lattice constants a = b = 7r plotted in the reduced zone scheme. 

In figure |21 we also show the sections of constant energy for the bands depicted in 
figure [T] By changing the amplitude of the potential the band structure will change, 
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and the bands become more flat and more separated as the amplitude of the potential 
is increased. In figures E0 we show the contour plots of Bloch states at different points 
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Figure 2. Sections of constant energy for the lower (a) and upper (b) bands in figure^ 
The white regions denote areas of higher energy. 

in the BZ. We remark that, due to the separability of the potential, both the derivative 
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Figure 3. Contour plots of the Bloch states at the corner k x = 1, k y = 1 of the BZ for 
the lower (a) and upper (b) band of figure 1. The wavefunction in (a) leads to soliton 
formation trough modulational instability while the one in (b) is modulationally stable. 



and the curvature of the band are independent on k y (respectively k x ) for fixed values of 
k x (respectively k y ) in the BZ. This is seen in figure El where the group velocity and the 
components of the reciprocal mass tensor are reported for the bands in figure ^ Similar 
calculations can be easily done for the 3D separable Mathieu potentials. The analysis 
can be extended to lattices with rectangular or tetragonal symmetry, as well as, to more 
general separable potentials such as the ones considered in |28j . 

In the following sections we shall use these results to compare numerical studies of 
the instability of Bloch states in presence of nonlinear interactions with the prediction 
(equation (JHJ)) of the previous section. 
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Figure 4. Contour plots of Bloch states of the lower band along lines of symmetry of 
the BZ. Case (a) refers to k x — 0.8, k y — 0.8 while (b) to k x = 1.0, k y = 0.8. 
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Figure 5. The group velocity (a) and the reciprocal of the effective mass uo xx 
(respectively w yy ) (b), as a function of k x (respectively k y ), for arbitrary k y 
(respectively k x ) in the BZ. The continuous and broken lines refer, correspondingly, to 
the first and second band of figure ^ (notice from (a) that there is no singularity for 
w XI (or Lo yy ) in the origin. 



4. Numerical simulations 

The linear stability analysis, described in section |21 gives the growth rates and spectra 
for the modulational instability. This information appears to be sufficient to predict the 
spatial arrangement and symmetry of the emerging soliton-like structures. However, to 
gain more insight into the development of modulational instability one has to recourse 
to numerical simulations. 

For the numerical study we have used the potential Vj(rj) = Acos(K,rj) for 
j = x,y,z, which is motivated by the recent experiments jHIH]- Then in the above 
formulas the terms corresponding to j = x, j = x,y, and j = x,y,z are retained, 
respectively for ID, 2D and 3D optical lattices. Also in the case at hand M"^, , 
and M~^, z have the same functional dependence on the arguments, which means that 
they coincide when q x = q y = q z . Numerical solution of the equation ((TJ has been 
performed by the operator splitting procedure using multi-dimensional fast Fourier 
transform j32J. The spatial domain x,y,z G [— §■••§] (i.e. L x = L y = L z = L) was 
represented by an array of 128 x 128 x 128 points. For ID and 2D cases the results 
were checked by increasing the number of grids (512 and 256 x 256 respectively), which 
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showed no qualitative difference. The time step was St = 0.001. To be specific, we 
concentrate on the case of positive scattering length, \ — 1-0, choose k = 2.0 (i.e. 
— CLy — &z — tt) , and p = 0.5. 

4-1. ID optical lattice 

Basic features of the development of modulational instability and formation of soliton- 
like excitations in effectively ID optical lattice was described in |%|ll2j. Below we extend 
the parameter values, which can lead to formation of qualitatively different types of 
localized excitations. 

The coefficient of nonlinearity x i n equation (JIJ is an important parameter which 
determines such a property of BEC as the macroscopic quantum self-trapping |2THlH0|IHlj . 
At strong nonlinearity the tunneling of atoms between adjacent wells of the optical 
lattice is suppressed, dispite the significant population imbalance, due to macroscopic 
quantum self-trapping effect. This property affects the development and further 
evolution of spatially localized excitations in BEC arrays. In figureElwe report two types 
of soliton-like excitations, developed in ID BEC arrays at weak and strong nonlinearities. 
All remaining parameters, except \i are similar in these two cases. Envelope soliton- 
like modes, which occupy few lattice sites are formed at weak nonlinearity, while the 
intrinsic localized modes, occupying a single lattice site are formed at strong nonlinearity. 
The time required for development of these excitations are also different. The localized 




Figure 6. Two types of soliton-like excitations in arrays of BECs confined to ID 
optical lattice, emerging due to the modulational instability of the initial waveform 
ip(x,0) = 0.5sin(x) at t = 0. The lower curve represents the optical lattice potential 
V(x) = Acos(kx). (a) Envelope solitons, occupying few lattice sites, are formed at 
t = 400 when the nonlinear atomic interactions are weak \ = 0.2. (b) Intrinsic 
localized modes, which fit a single lattice site, are formed at t = 70 for strong 
nonlinearity \ — 1-0- Parameter values A = 1.0, k — 2.0, L = 32tt. The initial 
waveform is perturbed by 6ip(x) = 0.01 sin(0.125x) 

excitations represented in the figure |B1 are thin disks filled in with BEC atoms, where the 
atomic density is much greater than in neighbouring array sites. The dynamics of these 
excitations is governed by the ID nonlinear Schrodinger equation jHj, and for particular 
parameter settings they can be stable ^2], or have very long (relative to duration of 
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experiments) recurrence times. The separability of the periodic lattice potential in 
equation (0) leads to similar scenarios of the development of modulational instability 
also in 2D and 3D cases considered below. 

4-2. 2D optical lattice 

Now let us consider in more detail the development of soliton-like excitations and the 
possibility to stabilize them in a 2D optical lattice. A 2D optical lattice is formed 
by overlapping two laser standing waves along the x and y axes, superimposed on a 
continuous BEC in a magnetic trap. The condensate is then fragmented and confined 
in many narrow tubes centered at lattice potential minima and directed along the z axis. 
As a result of modulational instability, the initial distribution of the atomic density over 
the tubes in the optical lattice is changed. 

In order to analyze the instability of initial waveforms we consider Bloch states 
corresponding to different points of the BZ. Let us consider the points qo = (±1,±1) 
at the boundary of the BZ (figure [TJ. Then, restricting consideration to the two lowest 
bands (a = 1,2), one can distinguish three different cases: 

4-2.1. Case 1. Both eigenfunctions \mo tX ) and \mo, y ) belong to the first lowest zone: 
m>o,x = m o, y = (1, ±1)- Then Mj"^ = = M^f 1 < and the wave is unstable. 

The BEC population dynamics in this case is reported in figure The most interesting 
feature of the modulational instability developed is that it evolves in a regular structure 
which represents symmetrically spaced localized in space (we call them soliton-like) 
distributions (see figure 0d)- Each of the humps shown in the figure represents a tightly 
confined tube along the z-direction. The number of tubes is proportional to the size of 
the box. In order to illustrate the last statement we performed calculations (see figure 
IHJ) with parameter settings similar to those of figure [7| with the exception of domain size 
L = 28tt. 

In order to understand this behavior we notice that from the equation (JHl) we get 
that the excitations with characteristic scales A > A m?r) = „ 27r = -\/ ^ M 3 - , where 

Tmn K m ax P V X ' 

M^ 1 is the inverse of the effective mass tensor, are unstable. The largest increment (i.e. 
the large Im|fi|) is achieved for A = y/2\ min . This has two consequences. First, the 
symmetry group of the developed structure must be of C n type with the symmetry axis 
coinciding with that of the condensate, and second, an effective scale A e // ~ Ao must 
be a characteristic scale of the most excitation which at the beginning of the evolution. 
One can estimate the value of Ao taking into account that for a chosen point of BZ the 
inverse of the effective mass tensor is (M]" 1 ! m 6 (see figure EJ) and for the solutions 
studied numerically the effective nonlinearity is x ~ 0.1935 (the respective normalized 
eigenfunctions are approximated by -sin(x) jH], which gives A ~ 20.053. This result 
corroborates with the distances between the humps along the radial direction measured 
from the direct numerical simulations: A e // ~ 23.0 (figureEJ). Next we have to take into 
account that the carrier wave mode is chosen at the point q = (±1, ±1) placed at the 




Figure 7. The evolution of BEC atomic distribution in a 2D optical lattice according 
to equation (JTJ with A = 1.0, K = 2.0, x = 1-0, L = 14n. (a) The initial waveform 
tp(x,y,0) — 0.5 sin(x) sin(y) at t = 0. (b) Formation of soliton-like excitations 
due to modulational instability at t=85 when the initial waveform is perturbed by 
5ip(x,y) — 0.01 sin(0. 125a;) sin(0.125y). (c) The distribution (b) slightly narrows 
but remains stable for very long times when the strength of the trap potential is 
adiabatically increased up to A = 3.0 during 85 < t < 90. The snapshot (c) is shown 
at t = 200 (stability is verified up to t=1000). 



corner of the BZ which corresponds to waves whose phases propagate in the directions 
x = ±y. This immediately specifies the symmetry C4. In other words, one can specify 
the points where the humps (confined tubes) should appear: in the plane (x, y) these 
are intersections of lines x = ±y with the circles of radii (| + pj A e // where p — 0,1, .... 
In a square box of size L one will observe L 2 /\ 2 e j^ humps. This estimate being rather 
rough (it does not take into account boundary effects) was confirmed by our numerical 
simulations. Also one can predict that the characteristic diameters of the humps should 
be less than X m in (~ 7.1 in our case). This gives an estimate for the BEC density in a 
tube n t versus the initial density n : n t = h^—n^, which in our case gives n t w 38no- 

min 

In order to evaluate the increase of the BEC atomic density in a soliton-like excitation, 
we have numerically integrated \ip(x, y, t)\ 2 at t = 85 (figure 0)) over individual lattice 
sites. The result is that 65 % of the BEC matter, initially uniformly distributed over the 
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Figure 8. Regular spatial structures with BEC in a 2D optical lattice, emerged due 
to modulational instability. Parameter settings are similar to those of figure except 
the domain size L — 287r. 



xxzxxn: 




Figure 9. Contour plot of localized excitations superimposed on the profile of the 
optical lattice potential, showing that each excitation fits a single cell. Dark regions 
correspond to the minima of the periodic potential for L x = L y = 12ir. 



optical lattice, are collected in four sites due to the modulational instability. Therefore, 
the increase of the atomic density in a localized excitation is rit = • 0.65 • uq ~ 32no, 
wich is close to the above analytical estimation. For the 3D case considered in the 
next subsection this last estimate for the BEC density in hollows rih is modified as: 

Lt, Lit L z 

n h = y z n . 

min 

As we have seen, the modulational instability results in formation of regular pattern 
of soliton-like excitations in arrays of BEC (figures [7|d, 8). However, they eventually 
decay in accordance with equation ©, which does not support stable solitonic solutions 
in 2D and 3D. A simple way to retain these excitations would be the increasing of the 
strength of the periodic trap potential, when excitations are formed. High potential 
barrier between lattice sites then suppresses the atomic tunneling, providing strong 
confinement. This idea is illustrated in figure where we show the evolution of the 
BEC atomic distribution. Until t = 85 the dynamics is guided by the modulational 
instability, resulting in formation of regular spatial structures with BEC. As the soliton- 
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Figure 10. Evolution of the BEC distribution corresponding to figure [7J shown as a 
diagonal section, d is a distance along the diagonal of the square domain with L = I4ir. 

like excitations are formed, the strength of the optical lattice was increased adiabatically, 
starting at t = 85 and ending at t = 90, according to law A(t) — 2[1 — | cos( 7r( - f ~2' f amp - > )], 
where t ram p = 85, At = 5. To obtain estimates in physical units recall that time is 
normalized to Jjp- = 50/xs (for 87 Rb atoms and ki = 8.06 ■ 10~ 6 m -1 ), and the strength 

h 2 k 2 

of the periodic potential is expressed in units of recoil energy E r = . Therefore, 
t = 85 corresponds to ~ 4 ms, while A = 1.0 corresponds to ~ E r . These values are 
typical in experimental situations. 

To understand the stabilization phenomenon we notice that X m in is of order of 
2a x {2a y ), which means that the most of BEC atoms are concentrated in a unique cell 
(see figure IHl). This type of excitations closely resemble the intrinsic localized modes in 
BEC in the tight-binding approximation [31H l34| . By increasing the potential amplitude 
one makes the optical lattice more deep, which in turn leads to decreasing both the 
probability of tunneling of atoms from the most populated cell to neighbor cells and a 
"number" of atoms in classically forbidden zone. As a consequence, the BEC density in 
the most populated cell is growing, which is illustrated by figures [?b and [TQJ 



4-2.2. Case 2. The eigenf unctions & m0yX and $ mi) !J belong to different zones, say $ mo x 
belongs to the first lowest zone: mo >x = (1, ±1) and $m , H belongs to the second lowest 
zone: m y = (2, ±1). Then < and > 0, and the condensate is unstable. 

In this case the instability condition takes the form < M^-K^ — iM^J-fTj < 

4%p 2 , and the most unstable excitations have K"i < i'f if 2 (which is related to the 

fact that an eigenfunction &m 0x belongs to the "unstable" branch). That is why the 
main instability results in a pattern having different symmetry: it develops in the re- 
direction. Along this direction the pattern is rapidly split in a sequence of solitary waves. 
The instability develops also along ^/-direction, but at much larger time scales (see figure 
ITTf . To estimate the number of humps, we take into account that the periodic boundary 
conditions impose a characteristic scale K x = 2ir/L which leads to the following estimate 
for the most unstable scale A , and thus to X e ff'- ~ 27r/ — |pT- For the case, 
illustrated in figure ITTk we obtain X pa 24, which yields for the number of humps in the 
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Figure 11. Distribution of BEC in a 2D optical lattice with A — 1.0, k = 2.0, x = 
1.0, L = 127r subject to the initial condition ■(/'(a;, y, 0) = 0.5 sin(cc) cos(y). (a) Regular 
pattern of soliton-like excitations are formed at t =32. (b) Contour plot of localized 
excitations showing that they occupy the sequence of single lattice sites along the y 
direction. 



x direction LK max /2n ~ L/A ~ 2. The instability is developed in ^-direction as well, 
which is characterized by much larger spatial and temporal scale, and for experimental 
purposes can be neglected. In order to preserve the developed spatial structures with 
high atomic density, it would be enough to increase the strength of the periodic trap 
potential, since the excitations fit the sequence of single lattice cells along the y direction 
(figure HUd). 

4-2.3. Case 3. Both eigenf unctions ®m 0x and &m 0iy belong to the second lowest zone: 
^0,2 = m o, y = (2,±1). Then M^.^ = M 2 ^ > and the wave is stable, which was 
confirmed numerically using the initial conditions ip(x,y,0) = 0.5cos(x) cos(y). 

4-2.4- Other points of the Brillouin zone. It is of particular interest for applications to 
explore the development of modulational instability of initial waveforms, corresponding 
to different points of BZ. We have tried the Bloch states with wavevectors spanning 
the first BZ. The main observation from the relevant numerical simulations is that, the 
modulational instability results in formation of different spatial structures with BEC 
depending on the initial waveform. The time required to formation of these structures 
is also different. 

As an example, in figure ^] we report the result of modulational instability of the 
Bloch state with mo tX = mo, y = (0.9, ±0.9). This initial distribution seems to be of 
interest because it resembles the situation, when a BEC of comparatively small size is 
loaded into the optical lattice, so that minor number of lattice sites are filled in with 
BEC, others being almost empty. In this case the cells in the central part of the optical 
lattice contain more BEC atoms, than the peripheral ones. The localized excitations 
developed due to the modulational instability occupy the central nearest-neighbor lattice 
cells. 
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Figure 12. Redistribution of the BEC atomic population in a 2D optical lattice due 
to modulational instability. A — 1.0, k = 2.0, \ — 1-0, L = Yin (a) The initial 
waveform, corresponding to the Bloch state mo^ x = rriQ ty = (0.9, ±0.9). (b) Formation 
of soliton-likc excitations due to modulational instability at t = 16. 



4-3. 3D optical lattice 

Qualitatively similar behaviour of the modulational instability with respect to formation 
of soliton-like excitations was observed in 3D case. The developed structures are small 
hollows filled in with BEC atoms of much greater density compared to surrounding 
array sites. Figure EH illustrates the emergence of spatial structures with high atomic 
concentration in a 3D BEC array, shown as a section along the main diagonal of the 
cubic domain with L = 12tt. The time interval is selected to display the emergence 
of soliton-like excitations at t ~ 28, and their subsequent decay. The relative atomic 




Figure 13. Formation of soliton-like excitations by t ~ 28 in a 3D BEC array, d is 
the distance along the main diagonal of the cubic domain with L = 12n. The initial 
condition is tp(x, y, z, 0) = 0.5 sin(a;) sin(y) sin(^). 



density in a localized excitation is estimated, as in the 2D case, by numerical integration 
of \ip(x,y, z,t)\ 2 at t = 28 over individual lattice sites. In the 3D case ~ 67% of all the 
BEC matter, initially uniformly distributed over 1728 lattice sites, are collected in eight 
sites due to the modulational instability, therefore rih = ■ 0.67 • n ~ 145no. 

For the considered initial parameter settings the long-term evolution of the atomic 
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distribution over the optical lattice in both 2D and 3D cases exhibits the recurrence 
phenomenon, which reveals itself as the return to the primary low density state. In 
order to prevent the decay of soliton-like excitations, similarly to 2D case, the strength 
of the periodic potential has to be increased adiabatically when the excitations are 
formed, e.g. at t ~ 28 for the above 3D parameter settings (figure US}. 

As far as the initial state is concerned, we remark that it could be experimentally 
realized starting from a uniform condensate with quasi-momentum k = (i.e. with 
equally filled potential wells), and accelerating the optical lattice in a particular 
direction. The acceleration of the lattice, being equivalent to a gravitational field, 
will induce the initial k = state to move along the band, until it reaches the edge 
of the band where the instability develops. Depending on the direction selected in the 
real space (this corresponding to a fixed direction in the reciprocal space), one will 
get different final localized states from the band edge instability. The issues relevant 
to preparation of particular Bloch states of the condensate in an optical lattice are 
discussed in p. 

It should be pointed out that the similarity of the development of modulational 
instability in all three dimensions is the result of separability of the periodic trap 
potential. Although this corresponds to the majority of experimantal situations, the case 
of non-separable potentials (non-orthogonal lattices) represents a significant interest. 

5. Conclusions 

We have studied the modulational instability in arrays of BEC confined to optical 
lattices. The formation of coherent spatial structures with BEC is shown to be the 
principal feature of the evolution of atomic distribution over the optical lattice, when 
guided by the modulational instability. In ID case the developed structures are the 
matter-wave solitons which may be regarded as thin disks of highly concentrated BEC 
atoms. Depending on significance of the nonlinearity, two distinct types of matter- wave 
solitons may develop, these being the envelope solitons (weak nonlinearity) occupying 
few lattice sites and intrinsic localized modes (strong nonlinearity), each fitting a single 
lattice site. In 2D and 3D cases the emerging spatial structures with BEC represent 
soliton-like excitations regularly arranged over the optical lattice which are, however, not 
stable. We proposed a simple way to stabilize these localized excitations by increasing 
the strength of the optical lattice when they are formed due to the modulational 
instability. Different initial waveforms, corresponding to particular points on the BZ 
are tested for instability. The aspects of developed spatial structures with BEC is 
shown to depend on the selected Bloch state. Theoretical model, based on the multiple 
scale expansion describes the primary features of emerging soliton-like structures with 
BEC, including the number of localized excitations, their spatial simmetry and relative 
density of BEC atoms they contain. The proposed method for creation and preserving 
of soliton-like spatial structures with highly concentrated BEC atoms, may be of interest 
for the physics and applications of BEC. 
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